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Abstract 

We present simulation results for single a-thermal chain polymers in finite volumes. For 
this we use a recently proposed recursive implementation of the enrichment method. In 
3 dimensions it allows the simulation of extremely long chains (up to N = 300,000). It 
is much less efficient for d = 2, but we can also there extend considerably the previously 
accessible range of chain lengths and densities. We verify all tested scaling laws except 
one, and we point out similarities with complex optimization problems. 
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1 Introduction 



Consider an ensemble of self avoiding chains ("SAW's") on a regular lattice in d 
dimensions. All chains start from the origin of the lattice. Each configuration of 
a chain of length N (here and in the following, N will be the number of links, so 
that the number of sites touched by the chain is iV+ 1) contributes to the partition 
function with a weight z , where z is a constant (iV-independent) "fugacity". If 
we denote by Cat the number of configurations with length N and with one end 
fixed at the origin (with Co = 1), we have thus 

oo 

Z(z) = £ z N C N . (1) 

N=0 

This provides a simple model for the grand canonical [] partition sum of an isolated 
polymer in a good solvent, where the excluded volume interaction is taken into 
account by the requirement of self avoidance. Although details are obviously not 
very physical, it should give the correct scaling behavior in the limit of very long 
chains, N — > oo. 

It is thus not surprising that an enormous amount of effort — both analytical and 
numerical |], — has been devoted to the study of the above model. In particular 
it is known that the sum in eq.(|l]) converges up to a critical value z c = 1/ fi, where 
Z{z) has a singularity 

„, . const , nS 

Z(z) ~ . 2 

K J (1 - nzp V ; 

Equivalently, Cn obeys a scaling law Cn ~ const /i^iV 7 " 1 . 

Most of the previous works have dealt with the sub critical phase z < z c , and the 
supercritical case z > z c has obviously been studied much less. In order to define 
the theory for z > z c , one has two alternatives: one can either make an analytic 
continuation, or one can regularize the sum in eq. ([I]) by a cut-off at large N. The 
latter can be done such that it has a very natural physical interpretation: we 
simply put the chain into a finite volume V, and count only those configurations 
which never leave the volume. 

In the limit V — > oo, we obtain thus a model for a single polymer chain which 
fills a macroscopic volume with finite monomer density p. In the literature this is 
known as the "dense limit" || |4]]. In the following we shall study this model by 
means of Monte Carlo simulations. 

Actually, the dense limit can also be understood as the limit of a more general 
model. There, the number of chains is fluctuating and controlled by a second 
fugacity x, and we take the limit x — > f|. There has been some controversy 

Notice that some authors H] call this the equilibrium ensemble, in order to distinguish it 
from a truly grand canonical ensemble in which also the number of chains can fluctuate; in the 
present case, we always work with a single chain. 
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in the literature following the claim by Gujurati || that there should be a phase 
transition in this model at a finite value of x. It seems now clear that this claim 
is wrong 0, ||, and we shall not mention this generalization further. In any case, 
the present paper will be devoted exclusively to the case x = 0. 

All relevant scaling laws at the transition from p = to p > can be obtained 
from the conventional (and well established) scaling theory for subcritical SAW's 
JT|. This will be reviewed in sec. 3, where also Monte Carlo data are given in 
support of these scaling laws in 3 dimensions. 

In the opposite limit, when p = 1 such that the SAW visits every point of the 
lattice, we are dealing with hamiltonian walks. These have been studied in several 
papers, the most comprehensive study being that of [0]. They form a model for a 
dense polymer melt. 

The intermediate region < p < 1 is particularly interesting in 2 dimensions, 
where it was shown in || that the model is a co-dimension ("self-organized" 
H) critical phenomenon. That means that various observables are described by 
non-trivial scaling laws, although there is no need for a control parameter to take 
any special value. In 3 dimensions there should be no such scaling laws, but we 
expect logarithmic behavior instead. 

It would of course be very desirable if one could verify these predictions numer- 
ically. Some numerical results based on transfer matrix calculations have been 
given m §. In that paper and in fll3| , |14| , also Monte Carlo simulations were 



presented. But only very small systems could be simulated in H], and systems 



with > 1 chain (i.e., x > in the above sense) were simulated in [13, [14]. Thus 



clear-cut and significant tests of the predictions of |3j would be most welcome. 
But all polymer Monte Carlo algorithms which we are aware of become inefficient 
at high densities. 

It is the aim of the present paper to show that recursive implementations || |TU| 
of the incomplete enumeration JXTJ] and of the enrichment |l!| methods can be 
quite efficient in simulating polymers in the dense limit. 

This is so in spite of the problem that we cannot justify our algorithms rigorously, 
in the way we have to apply them in most dense limit problems. While this seems 
to be mainly a formal problem in 3d (we do not encounter any serious problems 
unless we go to high densities), the situation is numerically very delicate in 2d. As 
a consequence we cannot go to really large systems and high densities, but we can 
go considerably much further than previous analyses. We will be able to check 
all scaling laws mentioned above with much higher accuracy than was possible 
previously. We will find good agreement in all cases except one. 

The main problem in 2d is that we cannot produce long chains if we from the 
beginning use a fugacity z which is much larger than the critical one. A better 
strategy is to start building chains with z very close to z c , and to increase it during 
the growth of the chain. Finding the optimal schedule which gives the longest 
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chains (i.e., the highest densities) is very similar to optimization in a problem 
with many local maxima, where greedy methods usually are sub-optimal. 

In sec. 2 we shall present our algorithm and the problems encountered when apply- 
ing it to dense systems. In sec. 3 we shall treat the limit of small densities, while 
our main results for scaling laws in 2c? systems will be discussed in sec. 4. The very 
dense limit and hamiltonian walks are treated in sec. 5. The paper concludes with 
a discussion in sec. 6. 



2 Algorithm 

We used a recursive algorithm which can be considered as a recursive and ran- 
domized implementation of the old enrichment method Jl2|, [K|]. In this algorithm 
we call a subroutine STEP(x,N) in order to add the next (iV-th) step to a chain 
starting at site x = and ending at x^r = x. Apart from updating whatever 
statistics we want to measure, this subroutine does the following: it first marks 
site x as occupied, calls itself a random number of times at one or more neigh- 
boring sites x' = x ± e, and with N replaced by N + 1 , and marks site x again as 
free before returning to the calling routine. The first call to the subroutine from 
the main program is with arguments (x , 1). 

In order to be concrete, we have to specify what we mean by "calls itself a random 
number of times" . Assume that the coordination number of the lattice is A/", so 
that we can continue a non-reversal walk without self avoidance at every step 
into M — 1 different directions. If we want to obtain a grand canonical ensemble 
(i.e. a sample in which each configuration is represented in average by const x z N 
chains), we have to make in average P = (A/" — l)z calls during each call of 
STEP, and each of the M — 1 neighbors must get in average the same number 
of calls. In incomplete enumeration JTT], [| this is achieved by taking independent 
probabilities for each possible continuation, i.e. for each free neighbor of x. In our 
present implementation, in contrast, we chose the neighbors such that the number 
of continuations shows the smallest fluctuation. More precisely, we first select a 
random non-reversal direction e and call STEP at the site x' = x + e, then we 
select randomly another of the non-reversal directions Q. If P < 2, we select a 
random number r G [0, 1] and make a second call iff r < P — 1. Otherwise, we 
choose a new random direction (different from all previous ones) and continue. If 
Q is the smallest integer < P, we try the first Q directions unconditionally, and 
use r only for the last trial. Obviously, the above assumes that P < J\f — 1. If 
P = M — 1, this algorithm is equivalent to exact enumeration. 

A conceptually trivial but technically very important generalization consists in 
letting P depend on N . If a sample obtained in this way contains nw,c chains 

2 In not too dense systems it is more efficient to skip the selection of a new direction, and to 
use the same random x' for all calls. 
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with configuration C, the contribution of C to the canonical partition sum is 



n N N 

As described here, enrichment is much more efficient than incomplete enumera- 
tion, typically by factors 5-30 depending on the lattice and on the density of chains. 
In order to understand why, we have to understand the temporal evolution of the 
chain length N. 

As pointed out in |D|, iV performs essentially a random walk with reflecting 
boundary at iV = 0. If z < z c , there is a bias for iV to decrease with (CPU- 
)time, since the average number of returns from subroutine calls is larger than the 
number of new calls. On the other hand, if z > z c , the walk is biased towards an 
increase of N. In either case, the algorithm will be inefficient in creating many 
statistically independent chains: in the first case since it does not create many 
chains at all, in the second because most chains will have common first parts and 
will thus be strongly correlated. The most efficient algorithm is obtained when 
z is such that there is no net bias. Notice that the latter is true also for dense 
systems. In this case, there is an additional bias towards smaller chains which 
increases with N and which is directly proportional to p. In an efficient algorithm 
we should thus start with P N « (Af — l)z c and increase P/v slowly with N in such 
a way that the net bias remains small. 

But absence of bias is not the only criterion for efficiency. The random walk 
in N is also characterized by an effective diffusion constant D. If D is small, 
then N changes little with time, and the algorithm is obviously inefficient. More 
precisely, the CPU time needed to create one new chain and erase it again (so 
that the next independent chain can be built up) is inversely proportional to D. 
For dilute systems one can easily convince oneself that D = 0(1) for incomplete 
enumeration, while it can be >> 1 for enrichment [I0[ |TB|, |T7|. The reason for 



the latter is easily understood. Assume that M » 1 so that there is very little 
attrition and we can choose P/v close to 1, more precisely P/v — 1 + 5. In this case 
most of the steps in the walk will be forward steps, but once we make a backward 
step (since we selected an occupied site for the next subroutine call) this step will 
have length 1/5 in average. Since the net bias is zero, we must thus have 1/5 
forward steps per each back steps, and D = 1/5. This is of course not rigorous 
(since the random walk is correlated), but it gives a roughly correct estimate. 

The main reason why our algorithms become inefficient at high density (even if 
P/v is chosen optimally, which is already not easy to achieve) is that the above 
argument break down there, and the effective value of D decreases dramatically. 
Again this is easily understood heuristically: in a very dense system, we have to 
make many attempts for each forward step, and the algorithm will spend much 
time in tempting to prolongate the chain in cases where there simply might be 
no possible prolongation. Notice that this cannot be avoided by adjusting P N so 
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that it would be more appropriate to each specific local configuration. has to 
be chosen uniquely for all local configurations, and has to be such that in average 
there is no bias. In this sense, polymers in the dense limit are even worse than 
SAW's in random media where we can chose a different optimized Pn(C) for each 
realization C of the randomness ||. 

In the following, we shall call "tour" the entire random walk between two returns 
to N = 0, i.e. the set of all chains generated between two successive returns to 
the main program. As described above, our algorithm has to be used such that 
statistics is collected while a tour is going on, but it must be read out only after 
a tour has been completed. Only in that case eq.(|3|) is correct. For instance, 
if we run our algorithm with P >> (A/ - — l)z c and stop it before finishing a 
tour, we will "override" some of the repulsive effect of self-avoidance constraint, 
and the end-to-end distance will in average be smaller than the correct one A 
crucial question for the following is whether it is possible (and if so, under what 
conditions) to disregard this problem, and to stop a simulation before a tour is 
finished. There are situations where this seems more or less the only feasible 
strategy, since otherwise we either get large finite size corrections (due to short 
chains which are included near the end points of a tour) or extremely large CPU 
times. 

Thus we would like to follow a strategy where we simulate essentially a single 
long tour at z > z c which would not be finished at any practically reachable time, 
and interrupt it after some prescribed number of subroutine calls. It seems that 
this works indeed surprisingly well in 3 dimensions (though we do not have any 
analytic proof for it), but it would fail dramatically in d = 2. 

The reason for the failure in d = 2 is easily understood. There, a typical chain 
gets stuck in traps as illustrated in fig. 1. If we call A the area of the trap, it takes 
~ e^ p ~ Pc "> A steps (with P c = (J\f — l)z c ) to get out of this trap again. Since the 
size of the trap can be of the order of the entire lattice, the logarithm of the time 
needed to get an unbiased sample is thus comparable to the logarithm of the time 
needed for a complete tour. 

In 3-d a chain does not form traps, and we expect this problem to be absent. 
Indeed, we found numerically that the chain length fluctuates during a tour only 
little around an average value which is independent of the random number se- 
quence (with a characteristic time of the density autocorrelation function which 
increases at most as a power of the system size), provided the number of steps 
is > V/(P — P c ). The latter is the number of steps which would be needed to 
generate a chain of length V on an infinite lattice. To illustrate the dramatic 
difference between the 2-d and 3-d cases, we show in fig. 2 histograms of chain 
lengths reached after 10 7 steps. Each histogram is based on > 100 independent 

3 The effective repulsion results from the fact that a tour which involves larger end-to-end 
distances will last longer and will thus contribute more to the average distance - provided we 
follow it until it is finished. 
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runs. In both cases, P was independent of N and 20% above the critical value, 
and the lattices had the same number of sites (512 x 512 resp. 64 3 ). While we see 
a very wide distribution centered at iV << V for d = 2 (panel a), the distribution 
for d = 3 (panel b) is sharply peaked at a finite density. This density is moreover 
found to be independent of the lattice size. 

The influence of trapping in d = 2 can also be seen in a different way. In fig. 3 we 
show the average chain lengths obtained after m = 10 k steps (k = 3, 4, 5, 6 and 7) 
as a function of z/z c (the average maximal length obtained during the first m steps 
shows very similar behavior, and is only slightly larger, as the chains hardly recede 
at z > z c \ we used z c = 0.3790524 [18[]). The most conspicuous feature in fig. 3 are 



the pronounced peaks at z/z c slightly above 1. The steep decay at z/z c < 1 is of 
course trivial. Less trivial is the fact that increasing z far above z c does not lead 
to longer chains. Instead it leads to more trapping. The algorithm becomes more 
"greedy" in the sense that there is more bias towards an increase of the chain 
length, and it cannot so easily leave a trap once it has run into it. In fig. 3 we 
show results from lattices of sizes up to 4096 x 4096 (which is practically infinite 
for 10 7 steps). In this limit of infinite lattices, the optimal value of z/z c seems 
to decrease towards 1 with the number m of steps. This means that greedyness 
is increasingly punished with chain length. The optimal value of z — z c scales 
roughly as 1/ y/m, and the chain lengths reached at these optimal values scale as 



The need for not being too greedy is very reminiscent of optimization problems 
involving complex functions with many local minima. Prototypes of such problems 
are provided by the traveling salesman problem, and by the ground state search 
in spin glasses ||19|| . The fugacity in the present case is analogous to the inverse 



temperature when treating such problems by simulated annealing |TJ|]. The main 
difference is that in our problem we do not have any frozen randomness. In typical 
hard optimization problems the multitude of false minima is due to randomness 
in external conditions which remain fixed during the optimization process. In the 
present case, in contrast, it is the initial stages of the process itself which were 
random and which remain "frozen" during a finite time which diverges however 
as the process goes on. In this respect our problem is closer to the complex 
optimization going on in evolving biological systems f2"D]| , where it is also mainly 
the complexity of the evolutionary process itself (the "co-evolution" ) which leads 
to multistability. 



3 Small Densities: Average Monomer Density 
versus Fugacity 
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3.1 Infinite Volume Limit 



The standard scaling theory for free SAW's |T], [!| describes the behavior of SAW's 
in the grand canonical ensemble in the limit z — > z c — e, i.e. if we approach the 
critical fugacity from the subcritical region. In the present section we propose 
that this scaling behavior can be extended slightly above z c . Thus we should be 
able to predict also the behavior of dense SAW's in the limit of small densities in 
terms of the usual critical exponents and scaling functions for non-dense SAW's. 

Let us consider the density a finite box of volume V = L d with periodic boundary 
conditions. Let p(x, z) be the density of monomers in an ensemble where a single 
polymer is attached with one end to the origin x = 0. The length of the monomer 
is determined implicitly by fixing the fugacity z, i.e. by giving each iV-step chain 
a weight oc z N . For 

e = inf£ < o (4) 

this density becomes independent of L for L —>■ oo, and obeys the scaling law |]| 

p(x, z) « r 1,v - d g{r 1/v e) , r = |x| , (5) 

where g(Q is a universal scaling function which can be assumed to be analytic 
at ( = 0. Notice that the latter is a direct consequence of our choice of scaling 
variable. More conventional would be the choice re u , but it is clear that this would 
generate a singularity in p(x, z) at z = z c if g were analytic. 

We now assume that eq.flU) holds also for e > 0. In this case the polymer should 
fill a large volume with constant density (except for a region near the origin), and 
thus p(x, z) —>■ const for x — > oo. This is only possible if 

9(0 ~ C^ 1 , (6) 

and from this follows |?T| 

Poo(z)= lim p(x,z) ~ e^ 1 . (7) 



In 3 dimensions we had no difficulty in checking this directly. We used systems of 
different sizes for different values of z, but for all z they were large enough so that 
the fluctuations in chain length were much smaller than the average chain lengths. 
This is necessary for the average density to be independent of the box size. Oth- 
erwise, we would have had a non-negligible chance that the chain length goes 
occasionally down to zero, which would obviously depend on box sizes. Specifi- 
cally, this meant for the runs closest to z c that we used box sizes 512 3 , and allowed 
chain lengths up to 3.3 x 10 5 . Notice that this implied also that we had to take 
statistics without waiting for tours to become completed, as discussed in the last 
section. 

Results from these runs are shown in fig. 4. In this plot we used z c = 0.2134908 
which is the value found in JT0] (with error ±0.0000005) for SAW's attached to a 
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planar surface. This value is slightly smaller than the previous best value by one of 
|22] for free SAW's in infinite volume, but additional unpublished simulations 



us 



with the latter geometry gave a value in perfect agreement with the above. In 
contrast, a much larger value with comparable accuracy (0.2134987 ± 0.0000010) 



was found from exact enumerations in |23||. We see from fig. 4 that we can again 



exclude the latter. Fig.4 is obviously compatible with scaling. A least square fit 
would give an exponent dv — 1 = 0.728 ± 0.1 or v — 0.579 ± 0.003, slightly below 
the best current estimate v = 0.587 |24|. Indeed, a careful look at fig.4 shows a 



slight downward curvature which might explain this small discrepancy. 

Before going on we should point out that eq.(|7]) has an important consequence 
for the free energy per monomer, as was pointed out in || . If we assume that the 
free energy is an extensive quantity in the thermodynamic limit, we should have 

±-lnC N , L ^-f(N/L d ) (8) 

with f(p) — > — In /i for p — > 0. If the grand canonical partition sum 

Z L {z) = £ z N C N , L « fdN zVe-WM (9) 

N J 

is estimated by the saddle point approximation, we find the saddle point at 

j-\pf(p)] = hz. (10) 
Putting f(p) = — ln/i + A/(p) = \iaz c + Af(p), we obtain 

-j- p [pAf(p)] = Mz c /z)~ p V(^). (11) 

Integrating this gives || 

f(p) = - ln/i + const pVC^-i) . ( 12 ) 



3.2 Finite Volume Corrections 

In d = 2 we cannot proceed in the same straightforward way for the reasons 
explained in sec. 2. Instead, we must work on lattices which are small enough, and 
with fugacities which are close enough to z c , so that all tours finish. Indeed, in 
order not to introduce any bias we have to prescribe the number of tours to be 
made, not the CPU time. 

This means however that we are never really in the limit described by eq.(^), and 
we have to consider finite-L corrections to it. Assuming that there is a unique 
diverging length scale £ ~ N u ~ e~ u in the polymer problem, any dependence 
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on L can only enter through the dimensionless ratio £/L, and we can make the 
ansatz 

p L (z) = (N)/L d « e^-VCeL 1 /") (13) 

where 0(C) is a new scaling function with 0(0) finite. Notice that this is only an 
ansatz for the spatially averaged density, and it implies no assumption about its 
x- dependence. 

Alternatively we can start with the canonical ensemble, i.e. with chains of fixed 
length N. Generalizing the usual ansatz for isolated chains in infinite volume and 
the discussion at the end of the last subsection, we assume 

C N , L « W-i e -"f(W d ) . (14) 

In contrast to eq.([L3|), the last ansatz should hold also for large densities, and 
even in the hamiltonian limit p = 1. In the latter, it is known that || 

c NiN * w^M" , (15) 

where 7# depends not only on the dimensionality d but also on the boundary 
conditions. In d — 2 it has been conjectured that 7# for all lattices is the same 
as for the Manhattan lattice, where it has the value 35/16 for periodic boundary 
conditions Q. In that reference it is also conjectured that all the power behavior 
in eq.(T5|) is for N < L 2 transformed into a term L 2<aH ~ l \ The latter would be in 



conflict with eq.(|Tj), whence we believe that it is just a trivial mistake. 



It is straightforward to show that eqs.(|13D and (|Tj) are mutually consistent for 
any value of 7 > 0, provided just that f(p) satisfies eq. (|T2"|) . Indeed, with 
eq.fllD the ansatz eq.(0) becomes C N)L « iV 7 " 1 /^ exp(-constiVp 1 /^- 1 )) = 
iV 7 " 1 p, x(N u I L) . Inserting this into the definition of (N) and replacing the 
summation over N by an integral gives eq. (|13"D . 



In order to verify eqs. (|bl| ) and (0) in 2 dimensions, we made simulations on square 
lattices of sizes between 16 x 16 and 256 x 256. During the simulations z was not 
kept fixed but depended on N (in order to reach longer chains), but as pointed out 
in sec. 2 it is trivial to compute Cn,l from this. In fig. 5 we plotted p(L) = e~^ p > 
as defined through eq.([b|). Actually, to reduce finite size corrections we replaced 
the argument of / by (iV — L)/L 2 . This is suggested by the fact that the finiteness 
of the lattice cannot be felt by chains with N < L (we used periodic boundary 
conditions). We see reasonable scaling, i.e. the functions / defined in this way 
really depend only on the chain density. We should stress that this would not be 
so if we had omitted the first factor in eq.([L4]). To verify eq. (|12"D we plotted Af 
in fig.6 against z — z c . We see that the prediction Af ~ p 2 is satisfied reasonably 
well. The agreement is certainly far from perfect, but it improves with increasing 
L, and it would be much worse if we would not have taken into account the first 
factor in eq.([L4|). 
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4 Spatial Properties in Dense 2-D Systems 



Some of the most striking predictions concern the spatial structure of dense SAW's 
in 2 dimensions. In spite of having no obvious control parameter which is set to 
any special value, the system is predicted to be critical, and several correlation 
functions should obey non-trivial anomalous scaling laws. We will not attempt 
here to rederive them but simply refer to for derivations. 



4.1 Grand Canonical End-to-End Distance Distribution 

The first saling law which we want to discuss here concerns the probability that 
a SAW has an end-to-end distance R, irrespective of the number of steps N it 
involves. Due to the effective repulsion provided by the self avoidance constraint, 
one should expect that this increases for small R. In the spin model analogy this 
means that the spin-spin correlation increases with distance, a paradoxical result 
which was put in question in . 



More precisely, we studied chains starting in the center of a circular disk of radius 
r cut out from a square lattice. No chain is allowed to reach further than r. For 



R << r, the prediction by Nienhuis [26| is 



P(R- z) ~ R 3/s . (16) 

In order to check this, we show in fig. 7 results from a disk with radius r = 45. 
The fugacity was chosen as z = 0.39, i.e. about 3% larger than the critical value. 
The data are normalized by dividing them by the number of chains with R = 1 
(which is also the number of self avoiding rings). We see in fig. 7 large fluctuations 
at low values of R which are not statistical, but due to the discreteness of the 
lattice. Statistical errors are much smaller. In the intermediate range of R we 
indeed observe the predicted scaling regime. The drop at very large R is simple 
due to the finiteness of the lattice (]. 

Analogous data for three dimensions are shown in fig. 8. There we show two curves, 
one for a moderately large spheric volume (r = 40) but reasonably large fugacity 
(z = 0.215 = 1.0071z c ), the other for a much larger volume (r = 77) but much 
closer to the critical point (z = .214 = 1.0023^ c ). Both curves show that P(R; z) 
increases for small R, but the increase is much slower than for d = 2 and does not 
seem to follow a power law. 



4 the choice of a circular boundary was mainly to keep this drop confined to a narrow region 
in R. 
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4.2 Number of Loops vs. Number of Open Chains 



Instead of working at fixed fugacity, let us now consider the number of SAW's 
and of self avoiding loops at fixed densities. Let us call the latter C^\- We work 
on square lattices with periodic boundary conditions, and count only loops which 
pass through the origin. According to ||, we expect 

r (o) 

~ L- 2 ^ , lD = 19/16 . (17) 

Numerically, this was checked in || for L < 16, at density N/L 2 = 0.4. Actually, 
in general 0.4L 2 is not integer. Thus for each L the chain length N nearest to 
0.4L 2 had to be taken, and the density slightly fluctuated around 0.4. While the 
global behavior was indeed as predicted, it was superimposed by strong odd/even 
oscillations. 

We found that this odd/even effect is mainly due to the following: due to the 
periodic boundary conditions our lattice is indeed a 2-d torus. Loops can either 
wind around the torus or have trivial topology (i.e., be homotopic to a point). 
The latter loops and those which wind twice around the torus always have N 
even. Walks which wind once around the torus have N even on lattices with even 
L, and vice versa. The total number of of loops depends smoothly on L. Thus the 
number of even-iV chains at even L are larger than the number of chains at L±l at 
any near-by value of N. In order to eliminate this effect, we thus plot in fig. 9 not 
precisely the ratio C^] L /C N ^ L but rather {Cf_ l>L + 2C { ff L + cP +1L )/2C N>L . We 
see indeed a perfectly straight line on a log-log plot up to L « 40 (the fluctuations 
for even larger L are statistical). The theoretical prediction is indicated in fig. 9 
by a straight line, and is obviously in perfect agreement with our data. 

Recently there has been some discussion about the value of 7^. Based on exact 
enumerations of collapsed self-attracting SAWs on unbounded lattices, Bennett- 
Wood et al. p7] found 7d = 0.92 ± 0.09, while exact solution for the honeycomb 



lattice with strip geometry gave 7d = 1 [28 . Our data clearly disfavor these 



values (our estimate is 7d = 1.14 ±0.05). One reason of the discrepancy with the 



work of [27] might be that collapsed self-attracting SAWs at finite temperatures 
are not equivalent to dense SAWs, in contrast to common belief |29| . Another 
reason might be that jo depends on the boundary condition, again in contrast to 
common belief. 



4.3 Number of Chains, Absolute Values 

As a last scaling law, we used the same data to check the prediction of || that the 
number of self avoiding chains at fixed density < p < 1 scales with a 7-exponent 
equal to 7^ + 1, i.e. as 

C>,l ~ L^e-^W (18) 
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(compare eq.([T5|). This would imply that the number of loops scales without any 
power prefactor. It agrees with the behavior in the hamiltonian limit, but we 
pointed already out in sec. 3 that it cannot be correct in the limit p — > 0. In the 
present subsection we shall test whether it holds for intermediate values of p, more 
precisely at p = 0.4. To get rid of the unknown constant prefactor in eq. (|i~8|) , we 
first form the ratio 

Cn(2L),2l/Cn(L),L ■ (19) 

Here N(L) is defined as the integer closest to 0.4L 2 . We then define an effective 
free energy per monomer as 

, f v _ ln[Cjv (2L)i2L /C JV(L)iL ] 

hip) ~ N(2L) - N(L) (20) 

and plot /l(p) against 1/L 2 . If eq.(|IB|) is correct, this should converge for L — > oo 
towards f(p) with slope [jp ln4]/[3p]. Otherwise, if eq.([L4]) is correct, the slope 
at 1/L 2 = should be [(7 — 1) ln4]/[3p]. This plot is shown in fig. 10. It suggests 
that eq. ([T8|) is not correct, and eq.(|I"4}) is preferred. A more precise statement 
about the value of 7 is difficult because of the large fluctuations visible in fig. 10. 
They are partly due to the fact that p is not strictly constant, but there seems 
also to be some large scale curvature superimposed. In view of this, our estimate 
from fig. 10 is 7 = 1.25 ± .10. This could mean that 7 indeed has its critical value 
43/32, and that eq.(O) holds for all finite densities. 



5 Very Dense Systems and Hamiltonian Limit 

It is clear from the above that our method becomes less and less efficient if we 
go to higher and higher densities. Thus it is clear that we can study hamiltonian 
chains only for very small systems. 

In fig. 11 we show results for 2-d (panel a) and 3-d (panel b) systems with free 
boundary conditions. Due to the boundary conditions, all statistics depend on 
the starting point of the chains. In fig. 11 we show average values, obtained by 
averaging over all L d starting points on a square (resp. cubic) lattice. The ob- 
servable shown is the decadic logarithm of the number of chains. Our results are 
in good agreement with expectations || |7j that 

hg(C N>L ) = \og[L- d C n,l] ~ -Nf(N/L d ) - B(N/L d )L d - 1 (21) 

starting pos. 

up to logarithmic terms. 

The lattice sizes shown in fig. 11 are rather modest: L < 10 for d = 2, and L < 5 
for d = 3. Indeed, we tried larger lattices, but we were simply unable to generate 
even a single hamiltonian chain on any larger lattice with our method. While we 
estimated that L = 6 for d = 3 and L — 11 for d = 2 should be feasible with 
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unduly large amounts of CPU time, lattices beyond that should be unaccessible 
with (variants of) our algorithm, and with present day computers. 

The data shown in fig. 11 were obtained by splitting the range [1, L d — 1] of chain 
lengths into 3 subintervals [l,Ni], [JVi + 1,JV 2 ] and [JV 2 + l,L d - 1]. On the 
middle interval we performed the Monte Carlo algorithm as described in sec. 2, 
while we made complete enumerations of all chains in the first and third intervals 
(typically, N x = 8 for d = 2 and N x = 6 for JV = 3, while N 2 « 0.75L d ). The 
exact enumeration on the first interval slightly reduced CPU time and statistical 
fluctuations, but was not very important for efficiency. In contrast, performing 
exact enumeration for N > N 2 was crucial. It implied that we found all extensions 
of random samples of chains with length N 2 . 

Indeed, as seen in fig.ll, the number of chains decreases dramatically when ap- 
proaching N — ► L d . For d = 2 and L = 10, e.g., we see that the number of 
hamiltonian chains is less than the number of chains with N « 75 by more than 
7 orders of magnitude! If we try to find hamiltonian chains by extending ran- 
domly chosen self avoiding chains which fill the lattice partially, we need thus 
> 10 7 chains of length 75 in order to have just a single success (the data shown 
in fig. 11a are based on 6 x 10 8 chains of length 75, those of fig. lib — where this 
factor is > 10 8 — are based on 1.4 x 10 9 chains). 

As seen from fig.ll, this decrease of Cjv,l when approaching the hamiltonian limit 
becomes rapidly stronger as we increase L, rendering our algorithm impracticable 
on larger lattices. Things are very similar for periodic boundary conditions. The 
situation is slightly better if we select the chains not randomly but such that dense 
chains are favored from the very beginning. Thus, the success rate is slightly higher 
if we use the Rosenbluth method , but not enough to treat lattices larger than 
11 x 11 resp. 6 



3 



We should stress that this represents a fundamental obstacle for any algorithm 
which attempts to generate hamiltonian walks by completing randomly chosen 
chains which partially fill the lattice. Just such an algorithm was presented in 
pnfl . In that paper first a non-random hamiltonian chain was formed (this is 



always easy). This chain is cut into two pieces at a randomly chosen site. The 
shorter half is erased, and completion of the longer part to a new hamiltonian 
chain is attempted using the Rosenbluth method. If this attempt is successful, 
the old chain is replaced, while it is restored in the case of failure. After that, a 
new random cutting site is chosen and the process is repeated. It was claimed in 



31f| that the algorithm was successful in generating random hamiltonian chains 
on 2-d lattices with up to 21 x 21 sites and with periodic be. According to our 
estimates this seems very unlikely. We conjecture that all successful attempts on 
the largest lattices in |TTJ involved cutting the chain near one of its ends. The 
chance for completing it when cutting near the center should have been practically 
zero, and thus the central piece of the chain should not have been changed at all. 

In order to improve the situation, one can use a sample of partial chains which is 



14 



strongly biased in favor of chains which do have at least one hamiltonian extension. 
We created such samples by growing random chains of length N rs L d /2 and 
discarding all chains for which the set of yet unvisited sites is disconnected. On 
square lattices with L — 11, e.g., this reduced a sample of chains with N = 60 
by a factor r* 10 3 . This rendered feasible complete enumeration of all extensions 
of the remaining chains. Repeating the test for connectivity of the unvisited sites 
several times improved further the efficiency In this way we could generate a few 
hamiltonian chains with L — 12 (d — 2) resp. L — 6 (d — 3), but on larger lattices 
the method failed again. 



Two surprising results were found in pi ): first of all, (R%) increased faster than 
linearly with N, (R%) ~ iV L26 , and secondly there were strong odd/even oscilla- 
tions superimposed on that growth. In fig. 12 we show our data for (R%), also 
for hamiltonian chains on square lattices with periodic boundary conditions. We 
see no indication for a faster than linear growth. We do see some odd/even os- 
cillations, but the amplitude is much weaker than that observed in pi] , and has 
opposite sign. 



6 Discussion 



In this paper we have presented simulations of single chains in the dense limit. 
We have seen that such simulations are very easy in > 3 dimensions when using 
a recursive and stochastic implementation of the enrichment method, unless one 
wants to go to very large densities. More interesting is the 2-d problem due do 
universal anomalous scaling laws which should hold there. We were indeed able 
to test these, with much higher precision than in previous analyses, although 
our method is much less efficient in 2d than in 3d. The reason for the latter is 
that our algorithm lets SAWs run into 'cages'. In the dense limit we have to 
use supercritical fugacity, implying that it becomes hard to leave a cage. This 
becomes the more pronounced the more the fugacity is above its critical value. 
This is very reminiscent of 'greedy' optimization algorithms in problems with 
frozen randomness. 

In spite of this problem we were able to perform simulations of much longer chains 
than previously. With one minor exception we verified all scaling laws predicted 
for dense 2-d polymers. This proves that this system is indeed a generic (co- 
dimension 1) critical phenomenon. 

We also performed simulations of the hamiltonian limit, and showed that neither 
our algorithm nor any variant thereof can be efficient. If one wants to create a 
random sample of hamiltonian chains by completing partially filling chains, one 
has to use for the latter a sample which is strongly biased in favor of such chains 
which have at least one hamiltonian extension. 

In the present paper we have reported only on simulations of a-thermal polymer 
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chains. In |T6[], we have simulated dense polymers near the #-point. These simu- 
lations were essential in eliminating surface effects in the coil-globule transition, 
and in understanding the nature of this transition. For more details we refer to 
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Figure Captions: 



Fig.l: Trapped 2-d SAW. In supercritical simulations using the present algorithm, 
there is a bias towards increase of the chain length, and it can take extremely long 
until such a trap is left again. 

Fig. 2: Histograms of chain lengths reached after 10 7 subroutine calls in super- 
critical simulations. Panel a: 2-d lattice of 512 x 512 sites, /i = 0.454863; panel 
b: 3-d lattice of size 64 3 , \i = 0.256189. 

Fig. 3: Average chain lengths obtained after 10 fc subroutine calls (k — 3, ... 7) on 

2- d lattices of sizes up to 4096 x 4096, plotted against z/z c . For each value of z, 
data are averages over 500 independent runs. 

Fig. 4: Log-log plot of the density of supercritical SAW's filling a finite but large 

3- d volume, in dependence of the fugacity /i. The critical fugacity was chosen 
according to JTD|] . 

Fig. 5: Effective coordination numbers per monomer for 2-d chains, after dividing 
out the logarithmic corrections from Cn,l (see eq. (|14"D . On the abscissa is plotted 
the corrected density (N — V)j L d . 

Fig. 6: Log- log plot of Af against the corrected density, obtained from the same 
data as in fig. 5. 

Fig. 7: Log-log plot of the number of chains with end-to-end distance R plotted 
against R 2 . These data were obtained on a disk of radius r = 45 at fixed fugacity 
z = 0.39. They are normalized to 1 at R = 1. The straight line has slope 3/16, 
and indicates the theoretical prediction. 

Fig. 8: Simular to fig. 7 but for three dimensions. 

Fig. 9: Number of loops divided by number of chains, both at fixed density p = 
0.4, plotted against L. The theoretical prediction is indicated by the straight line 
with slope -19/8. 

Fig. 10: Free energy per monomer fi(p = 0.4) against 1/L 2 . Indeed, since 0.4L 2 is 
generally not an integer, the density slightly oscillates around the value 0.4. This 
is mainly responsible for the fluctuations seen in the data, as statistical errors are 
comparable to the size of the dots. 

Fig. 11: Decadic logarithms of the number of chains on lattices with free bound- 
aries, averaged over the position of the starting point. For d = 2 (panel a), curves 
are for L = 4, . . . 10, while they correspond to L = 3, 4 and 5 for d = 3 (panel b). 

Fig. 12: (i? 2 ) 1 / 2 against L for hamiltonian chains on square lattices with periodic 
boundary conditions. The data for L < 6 are based on exact enumerations and 
thus free of statistical errors. 
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